# Define priors for Baysian data analysis and conduct prior predictive simulations

# Load necessary libraries
```{r}
library(rethinking)
```

# Load dataset
```{r}
# Load data from module_1
load("./data/d_sol_cw_w_in.RData")

# Load simulation data from module_2
load("./data/d_sim_in.RData")

d_in <- d_sol_cw
#d_in <- d_sim_m1
#d_in <- d_sim_m2
#d_in <- d_sim_m3
```

# 1. Find prior for Sol -> T_S
```{r}
n <- 100
a_T <- rnorm(n, 0, 1)
b_T_S <- rnorm(n, 0.5, 0.5)

windows();
plot(NULL, xlim = range(d_in$Sol_std), ylim = range(d_in$T_S_std),
xlab = "Sol (std) [W/m^2]",
ylab = "T_S (std) [K]")
xbar <- mean(d_in$Sol_std)
for (i in 1:n) {
    mu <- a_T[i] + b_T_S[i] * (d_in$Sol_std - xbar)
    lines(d_in$Sol_std, mu, lwd = 1, col = 4)
}
points(d_in$Sol_std,d_in$T_S_std,col="black", lwd= 1)
```

# 2. Find prior for Sol -> C_MW
```{r}
n <- 100
a_C <- rnorm(n, 2.39, 1)
b_C_S <- rnorm(n, 0.8, 0.5)

windows();
plot(NULL, xlim = range(d_in$Sol_std), ylim = range(d_in$C_MW),
xlab = "Sol (std) [W/m^2]",
ylab = "C_MW [pg/L]")
xbar <- mean(d_in$Sol_std)
for (i in 1:n) {
    mu <- a_C[i] + b_C_S[i] * (d_in$Sol_std - xbar)
    lines(d_in$Sol_std, mu, lwd = 1, col = 3)
}
points(d_in$Sol_std,d_in$C_MW,col="black", lwd= 1)
```

# 3. Find prior for T_S -> C_MW
```{r}
n <- 100
a_C <- rnorm(n, 2.39, 1)
b_C_T <- rnorm(n, 1, 0.5)

windows();
plot(NULL, xlim = range(d_in$T_S_std), ylim = range(d_in$C_MW),
xlab = "T_S (std) [K]",
ylab = "C_MW [pg/L]")
xbar <- mean(d_in$T_S_std)
for (i in 1:n) {
    mu <- a_C[i] + b_C_T[i] * (d_in$T_S_std - xbar)
    lines(d_in$T_S_std, mu, lwd = 1, col = 2)
}
points(d_in$T_S_std,d_in$C_MW,col="black", lwd= 1)
```

# 4. Find prior for W -> C_MW and W -> T_S
```{r}
n <- 100
a_C <- rnorm(n, 2.39, 1)
a_W <- rnorm(n, 0, 1)
b_C_W <- rnorm(n, 0.5, 0.5)
b_T_W <- rnorm(n, 0.5, 0.5)

windows();
plot(NULL, xlim = range(d_in$W_std), ylim = range(d_in$C_MW),
xlab = "W (std) [m/s]",
ylab = "C_MW [pg/L]")
xbar <- mean(d_in$W_std)
for (i in 1:n) {
    mu <- a_C[i] + b_C_W[i] * (d_in$W_std - xbar)
    lines(d_in$W_std, mu, lwd = 1, col = 3)
}
points(d_in$W_std,d_in$C_MW,col="black", lwd= 1)

windows();
plot(NULL, xlim = range(d_in$W_std), ylim = range(d_in$T_S_std),
xlab = "W (std) [m/s]",
ylab = "T_S [K]")
xbar <- mean(d_in$W_std)
for (i in 1:n) {
    mu <- a_W[i] + b_T_W[i] * (d_in$W_std - xbar)
    lines(d_in$W_std, mu, lwd = 1, col = 3)
}
points(d_in$W_std,d_in$T_S_std,col="black", lwd= 1)
```

# 5. Find prior for Sol -> r_W, T_S -> r_W, r_W -> W, and r_W -> C_MW
```{r}
n <- 100
a_C <- rnorm(n, 2.39, 1)
a_r <- rnorm(n, 0, 1)
b_r_S <- rnorm(n, 0.5, 0.5)
b_r_T <- rnorm(n, 0.5, 0.5)
b_r_W <- rnorm(n, 0.5, 0.5)
b_C_r <- rnorm(n, 0.5, 0.5)

windows();
plot(NULL, xlim = range(d_in$r_W_std), ylim = range(d_in$C_MW),
xlab = "r_W (std) [l/min]",
ylab = "C_MW [pg/L]")
xbar <- mean(d_in$r_W_std)
for (i in 1:n) {
    mu <- a_C[i] + b_C_r[i] * (d_in$r_W_std - xbar)
    lines(d_in$r_W_std, mu, lwd = 1, col = 3)
}
points(d_in$r_W_std,d_in$C_MW,col="black", lwd= 1)

windows();
plot(NULL, xlim = range(d_in$r_W_std), ylim = range(d_in$T_S_std),
xlab = "r_W (std) [l/min]",
ylab = "T_S [K]]")
xbar <- mean(d_in$r_W_std)
for (i in 1:n) {
    mu <- a_r[i] + b_r_T[i] * (d_in$T_S_std - xbar)
    lines(d_in$T_S_std, mu, lwd = 1, col = 3)
}
points(d_in$T_S_std,d_in$r_W_std,col="black", lwd= 1)

windows();
plot(NULL, xlim = range(d_in$r_W_std), ylim = range(d_in$Sol_std),
xlab = "r_W (std) [l/min]",
ylab = "Sol (std) [W/m^2]")
xbar <- mean(d_in$r_W_std)
for (i in 1:n) {
    mu <- a_r[i] + b_r_S[i] * (d_in$Sol_std - xbar)
    lines(d_in$Sol_std, mu, lwd = 1, col = 3)
}
points(d_in$Sol_std,d_in$r_W_std,col="black", lwd= 1)

windows();
plot(NULL, xlim = range(d_in$r_W_std), ylim = range(d_in$W_std),
xlab = "r_W (std) [l/min]",
ylab = "W (std) [m/s]")
xbar <- mean(d_in$r_W_std)
for (i in 1:n) {
    mu <- a_r[i] + b_r_W[i] * (d_in$W_std - xbar)
    lines(d_in$W_std, mu, lwd = 1, col = 3)
}
points(d_in$W_std,d_in$r_W_std,col="black", lwd= 1)
```

## Run all above
```{r}

```
